

saux = (((1-alpha)*s*lambda)/((sigma^2)/2))  -  (1/(1-alpha));
Ahi  = (saux/(mL^saux)) / ( ((mM/mL)^saux) - 1 );
Alo  = Ahi/100;
obs_const= 100; 
constvec = linspace(Alo, Ahi, obs_const)';   % 1 x obs_const


%% Natural wastage
Gfvec_mM = (constvec/saux)*(mL^saux)*( (mM/mL)^saux - 1 );
gfvec_mM = constvec*(mM^(saux-1));


%% Hiring
neintegral = zeros(obs_ne,1);
for i = 1:obs_ne 
    if i >1
        nu = mu - (1-alpha)*( eta_ne(1:i) - quit_ne(1:i) );
        neintegral(i) = wghts_ne(1:i)'*((nu - sigma^2)./grid_ne(1:i));
    end
end
gfvec_ne = exp( (2/(sigma^2))*neintegral );
Gfvec_mU = wghts_ne' * gfvec_ne;
Gfvec_mU = Gfvec_mU*gfvec_mM; 


%% Solve for unknown:
Integral = Gfvec_mM + Gfvec_mU;
constsol = interp1(Integral, constvec, 1);


%% Final distribution function
gf_mM = constsol*(mM^(saux-1));
Gf_mM = (constsol/saux)*(mL^saux)*( (mM/mL).^saux - 1 );
Gf_nw = (constsol/saux)*(mL^saux)*( (grid_nw/mL).^saux - 1 );
Gf_ne = Gf_mM + gf_mM*cumsum( wghts_ne .* gfvec_ne);
Gf_all = [Gf_nw; Gf_ne];


